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Abstract 

We study, using noise reduction techniques, layer by layer epitaxial growth in 
limited mobility solid-on-solid nonequilibrium surface growth models, which 
have been introduced in the context of kinetic surface roughening in ideal 
molecular beam epitaxy. Multiple hit noise reduction and long surface dif- 
fusion length lead to qualitatively similar layer by layer epitaxy in 1+1 and 
2+1 dimensional limited mobility growth simulations. We discuss the dy- 
namic scaling characteristics connecting the transient layer by layer growth 

regime with the asymptotic kinetically rough growth regime. 
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I. INTRODUCTION 



Thin film growth, under solid-on-solid epitaxial conditions, from the vacuum vapor de- 
position of an atomic or a molecular beam (the so-called molecular beam epitaxy or MBE) is 
an important technological process used extensively to produce high quality thin films with 
smooth and flat surfaces and interfaces. It is also a growth process of considerable funda- 
mental significance [|T]-0] in the statistical mechanics of nonequilibrium phenomena because 
MBE (at least in its ideal form |7|], with no evaporation and vacancy or overhang formation 
at the growth front) in principle represents U\M a novel universality class of nonlinear surface 
growth outside the generic Kardar-Parisi- Zhang universality ||. A great deal of attention 
has therefore focused over the last ten years on the statistical properties of kinetically rough 
(and, in principle, generically scale invariant) surface growth in low temperature (room tem- 
perature or below) MBE, following the suggestions |7|,|8|,p!0|- [l^1 of the possible importance of 
ideal MBE in defining novel growth universality classes in kinetic surface roughening |l] Q . 
We note that the conserved surface current nature ]7||| of ideal MBE growth (i.e. solid- 
on-solid growth with no evaporation and vacancy or overhang formation) rules out a KPZ 
description of its growth dynamics. 

It is interesting, perhaps even ironic, that MBE growth has played such a central role 
]T|-^|] in kinetic surface roughening phenomena because the primary materials science impe- 
tus for MBE growth is obviously to avoid kinetic roughening as completely as possible in 
order to produce smooth and flat thin films of high surface quality with minimal amount of 
surface roughness. From the materials science perspective of producing high quality smooth 
(i.e. manifestly non-rough) thin films, therefore, MBE is typically carried out at elevated 
temperatures (~ 500-1000 K), where fast surface diffusion enables one to produce smooth 
thin films with very little surface roughness. Smooth MBE growth, as opposed to kineti- 
cally rough (low temperature) growth, is characterized by layer-by-layer growth oscillations 
Hl3|Jl4}| , where each layer of the growing thin film (on a singular high symmetry substrate) 
essentially fills up completely before the next layer deposition begins (on the other hand, 
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kinetically rough growth is, by definition, multilayer as many layers are partially filled at 
the growth front producing increasing surface roughness with no layer-by-layer oscillations), 
consequently the surface morphology and associated properties oscillate as growth proceeds. 
In the ideal layer-by-layer growth mode, therefore, the interface width (W, the root mean 
square fluctuation in the interface height) of the growing film oscillates (nominally between 
and 1, as measured in lattice units, indicating an empty or a filled layer) whereas in the 
kinetically rough growth mode W increases monotonically as a power law in the average 
film thickness {(h)). Layer- by-layer epitaxy is, however, an initial transient growth regime 
which eventually crosses over to the asymptotic kinetically rough growth regime as the shot 
noise intrinsic in the incident deposition beam fluctuations always wins out to damp the 
layer-by- layer oscillations, and at long enough time scales (and for large enough lateral sys- 
tem sizes) statistically scale invariant kinetically rough growth would always emerge. (In 
fact, the noise associated with the stochastic diffusion process also contributes to kinetic 
roughening, but the shot noise associated with the incident beam fluctuations is the most 
relevant roughening mechanism.) This is also the experimental observation : layer- by-layer 
growth oscillations, as studied for example through RHEED intensity oscillations monitoring 
the dynamical surface evolution [|D3|,[rjJ], eventually always damp out as the stochastic de- 
positopn shot noise associated with incident particle beam fluctuations leads to kinetically 
rough multilayer growth after some characteristic time t c . The damping time t c , beyond 
which layer by layer growth dies out, depends on the growth temperature (which controls 
the surface diffusion rate), and is in general larger for higher temperatures because longer 
diffusion lengths at higher temperatures enhance layer by layer growth. (There is actually 
a complication, arising from the unaviodable vicinality in the starting substrate which can 
never really be precisely a high symmetry singular plane in real growth where layer by layer 
oscillations tend to disappear at both high and low growth temperatures — the low tem- 
perature behavior is from the multilayer kinetically rough growth as discussed above, but 
the high temperature disappearance arises from the so-called step flow growth mode which 
is caused by the very fast surface diffusion of deposited atoms at high temperatures leading 
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to their moving directly to step edges, which must be present in any real substrate due to 
vicinality or miscut, without any layer by layer growth oscillations; we neglect considerations 
of such a step flow growth mode in this paper assuming all growth to be occuring strictly 
on singular high symmetry substrates.) 

The ideal MBE growth (on flat singular high symmetry substrates) can be thought of as 
composed of two regimes — an early time (t < t c ) transient (fast diffusion driven) regime 
of layer by layer growth followed by the asymptotic (t > t c ) kinetically rough (deposition 
beam fluctuation driven) growth regime (with no layer by layer oscillations) characterized 
by power law evolution [|IHL2|1 in surface roughness. At low enough temperatures, when 
surface diffusion is extremely slow, t c could be less than the time it takes to grow one 
monolayer of deposit on the average, and in that situation layer by layer transient growth 
regime is invisible with the kinetically rough growth regime being dominant essentially right 
from the beginning. Conversely, sufficiently high temperature growth on a small substrate 
could continue in the layer by layer mode for a very long time although some damping of 
the growth oscillations is inevitable with time as distant spatial regions on the substrate 
must lose coherence due to the inherent shot noise fluctuations associated with the discrete 
deposition process in the incident beam. Thus layer by layer growth us purely a "finite size" 
(both spatially and temporally) transient phenomenon — if the substrate is made sufficiently 
large and/or if one waits for sufficiently long growth time, layer by layer epitaxy must 
necessarily cross over to kinetically rough growth. It is important to emphasize, however, 
that the surface diffusion length is typically an exponentially activated function of growth 
temperature, and therefore a small change in temperature could cause a sharp large change 
in the growth mode (layer by layer to rough and vice versa depending on whether the 
growth temperature is decreased or increased) for a given substrate, leading to the empirical 
concept [[HJ of an epitaxial growth temperature T c with growth being layer by layer (rough) 
for T > T c (T < T c ) — clearly T c is a loosely defined concept because it must be a weak 
(sub- logarithmic) function of the effective substrate (or the terrace) size for a given material 
|15| . In general, "good" MBE growth aiming toward producing high quality smooth epitaxial 
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thin Aims is carried out at the highest possible growth temperature (within the constraint 
that evaporation or desorption from the growth front should be negligible so that the growth 
temperature cannot be arbitrarily high) so as to make atomic mobility at the growth front to 
be very high leading to large "surface diffusion length" I. Here I is taken to be the linear size 
over which the surface is smooth due to atomic diffusion. Assuming the deposition process 
to be a random Poisson process it is then easy to see that the typical surface roughness over 
terraces of size / would only grow as J (h)/l, where (h) is the thickness of the grown film (and 
we measure all lengths in lattice units). Thus for large /, one would have to grow a very thick 
film of thickness I 2 before the surface roughness reaches even one monolayer fluctuation. One 
can, therefore, grow MBE thin films of very high smoothness and quality, without worrying 
at all about the kinetic surface roughening by properly adjusting the growth temperature 
|L5] to make / large. 

Layer by layer MBE growth has been extensively studied ||T6|-[20|] in the literature us- 
ing computer simulations of MBE growth through the stochastic (or kinetic) Monte Carlo 
simulations, with the atomistic diffusion at the growth front assumed to be controlled by 
stochastic activated hopping process with the hopping rate determined by a local coordi- 
nation number dependent activated Arrhenius hopping. Such activated diffusion Arrhenius 
hopping simulations (sometimes also referred to as "full diffusion" simulations to differenti- 
ate them from "limited mobility" growth models which are our main interest in this paper) 
involve continuous (possible) hopping of all surface atoms according to their local bonding 
configurations (which determine the activation energy for the hopping process). Such full 
diffusion simulations are obviously not well designed to study the kinetic roughening univer- 
sality class of MBE growth because they are extremely time consuming and cannot really 
be carried out for large systems (particularly in the physically relevant 2+1 dimensions) 
for long times, an essential requirement for ascertaining the asymptotic universality class 
of a growth model. Although there are some notable exceptions [f2~T1 -|2"3|j , the full diffusion 



Arrhenius activated kinetic Monte Carlo simulations of MBE growth have not been used 
with particular success for understanding statistical scale invariance properties of kinetic 



surface roughening. Instead, important insights into the MBE universality class of kinetic 
surface roughening have come primarily from nonequilibrium limited mobility growth mod- 



els — mainly the so-called Wolf- Villain |]TT| (WV) and the Das Sarma-Tamborenea [10 
(DT) model — which were introduced specifically for the elucidation of the MBE growth 
universality. 

In this paper we study the DT and the WV model (we emphasize that WV and DT 
models, in spite of their close similarity in growth rules, belong |24[] to different asymptotic 
universality classes in both 1+1 and 2+1 dimensions although their pre-asymptotic scaling 
behavior is very similar which has led to considerable confusion in the literature) in the 
complementary layer by layer growth regime rather than the kinetically rough growth regime 
which motivated the introduction of these models. We mention in this context that some 
1+1 dimensional studies of WV and DT models in the layer by layer growth regime have 



recently been reported in the literature [p5| , p6|| . Our results, where applicable, agree with 



these earlier works |25|,|26|], but our focus in this paper is 2+1 dimensional growth and the 
effect of long surface diffusion length in 1+1 dimensional growth, neither one of which has 
earlier been studied. 

In limited mobility growth models (the models and growth rules used in this paper are 
described in section II of the paper — see, for example, Fig. []]), in sharp contrast to full diffu- 
sion MBE growth simulations, the goal is to suppress crossover and transient effects as much 
as possible (so as to efficiently reach the asymptotic kinetic surface roughening regime) and 
as such only the most recently deposited atom is allowed to diffuse or relax instantaneously 
to the appropriate incorporation site following the mobility rules of the specific model. This 
allows suppression of crossover effects invariably present in the full diffusion simulations 
arising from many different diffusion rates corresponding to many different possibilities for 
local bonding configurations — the only time scale in the limited mobility growth models 
being the deposition rate, which defines the time unit for the problem. From now on we take 
the time unit (sometimes referred to as a "second" ) as the time to deposit one monolayer on 
the average. Thus the growth time in this paper also defines (h) — the average thickness of 
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the deposited film measured in units of monolayers or the lattice constant, which we take to 
be the unit of length throughout. (With no loss of generality we choose the lattice constant 
to be the same along the substrate and the growth directions.) 

The limited mobility growth models [rti|,|n| are by construction strongly dominated by 
the deposition shot noise because the goal is to study the scale invariant kinetic surface 
roughening behavior. This is particularly true in the original versions of the growth model 
where the surface diffusion length is choosen to be unity, I = 1, i.e. the deposited atoms 
are allowed to move only to the nearest neighbor incorporation sites around the deposition 
site. The original DT and WV models therefore did not exhibit, by design, any layer by 
layer growth oscillations since the smoothing or the healing distance (7) was just one lattice 
unit. In order to manifest layer by layer epitaxy in limited mobility growth models one must 
therefore suppress the shot noise associated with the incident beam fluctuations. 

In this paper we accomplish the noise suppression by two alternative techniques : The 
'multiple hit' noise reduction technique [p7| -|3T[] and the 'long surface diffusion length (/ > 1)' 
noise reduction technique J2TJ. These techniques, described in section II of the paper, give 
rise to layer by layer growth (as monitored by an oscillatory surface roughness, i.e. W(t) 
showing oscillations as a function of growth time t) in the limited mobility growth models 
as described in section III of the paper. 

The rest of this paper is organized as follows. In section II we describe the limited 
mobility growth models and the noise reduction technique (s) employed by us. We also 
provide some theoretical background for our analysis of the simulation results. In section 

III we present and describe our numerical simulation results for layer by layer epitaxial 
growth in DT and WV models. We also discuss in section III various (approximate) scaling 
properties of our simulated layer by layer epitaxial growth results. We conclude in section 

IV with a general discussion of our results making connections with some of the existing 
results in the literature and pointing out possible future directions as well. 
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II. MODELS, THEORY, AND BACKGROUND 



The DT [1C] and WV [1 1| models used in our simulations are shown in Fig. |I[ We carry 
out growth simulations in both 1+1 dimensions and 2+1 dimensions (on 100 high symmetry 
substrates). Particles are dropped, obeying solid on solid constraint, sequentially (one by 
one) and randomly with an average rate of N per second, where N = L for d=l+l and 
N = L 2 for d=2+l, on a substrate of lateral linear size of L lattice units. We measure length 
in lattice units and time in inverse deposition rate (i.e. in units of monolayer filling since one 
average monolayer of deposition occures every "second"). Each deposited atom is allowed to 
"diffuse" instantaneously to its incorporation site following the mobility rules of the specific 
model. The diffusion rules in the DT model are that a deposited atom can move only if it 
has no lateral nearest neighbor in the same layer (if it does, then the atom is incorporated at 
the deposition site) — if the deposition site has no lateral nearest neighbor then the incident 
atom may move instantaneously to a neighboring empty site (within a lateral diffusion length 
of I, where I = 1 in the original model and in most existing simulations) provided the final 
incorporation site has a higher lateral coordination number (i.e. one or higher) than the 
deposition site. If several neighboring sites satisfy the diffusion rule then the atom will 
move randomly to any one of them with equal probability. The rules for the WV model are 
superficially similar to the DT rules : In the WV model all deposited atoms (and not just 
the ones with no lateral bonds) can, in principle, move provided they can increase their local 
lateral bonding and the deposited atom always moves to the site with the maximum local 
bonding environment. In both models, the deposited atom is incorporated at the deposition 
site if it cannot satisfy the diffusion rules (i.e. no sites with higher coordination available) 
within the diffusion length. 

Both DT and WV models have been extensively studied in the literature (mostly within 
nearest neighbor, 1 = 1, diffusion rules) in the context of their kinetic surface roughening 



universality classes. Recently, layer by layer epitaxy in the WV p5| , |26| , |30| and the DT 
PH model have been investigated in 1+1 dimensions using the multiple hit noise reduction 
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technique. The very first simulational observation of layer by layer growth in a limited 
mobility growth model was reported in the DT model in Ref. |21[ where it was studied in 
1+1 dimensions using a long (7 > 1) surface diffusion length, but no details were investigated. 
We emphasize that the usual I = 1 limited mobility growth model does not exhibit any layer 
by layer epitaxy by definition, and manifests kinetic roughening right from the beginning 
since for I = 1 the layer by layer epitaxy regime is restricted to less than one monolayer 



coverage, i.e. in the standard limited mobility growth models [fUJlIU the layer by layer 
epitaxy regime does not exist. 

To obtain a layer by layer growth regime in the DT and the WV model we use two distinct 
techniques to suppress noise and enhance diffusion, which enable our growth simulations to 
manifest strong layer by layer growth oscillations before crossing over to the kinetic surface 
roughening regime with pure multilayer growth. One technique, referred to as the noise 
reduction technique, has also been used by others p5|-p0|] to produce layer by layer growth 
in various models, e.g. Eden model ||28|| , single step model |29[, and the WV and the DT 



model in 1+1 dimensions [F2^,^B|,[J(J . We have earlier used this technique J3TJ to suppress 
corrections to scaling in the asymptotic kinetically rough surface growth regime in the DT 
and the WV model in order to accurately determine the dynamic scaling exponents and 
the associated growth universality class. In the noise reduction technique, characterized by 
an integer number m (the usual growth model without any noise reduction is an m = 1 
model), a counter is put on each surface site and each discrete deposition event on a site 
advances the counter by unity. A deposition event at a particular site is accepted only when 
the counter reaches a pre-determined number m > 1. Thus, this technique is the multiple 
hit noise reduction technique since m(> 1) deposition hits on a site are needed for a true 
deposition. After m hits on a site (i.e. after the acceptance of a deposition event) the counter 
at that particular site is set back to zero, and the whole multiple hit process begins all over 
again. The multiple hit noise reduction technique is a coarse-graining procedure which 
suppresses the deposition shot noise, and the noise reduction is enhanced for larger values 
of m. The second technique applied by us for obtaining layer by layer growth oscillations in 



limited mobility growth models is to use long surface diffusion lengths (/ > 1) in the growth 
simulations. Obviously long diffusion lengths enhance the layer by layer growth regime, and 
in particular, for I > L i.e. the surface diffusion length exceeding the system size, the layer 
by layer growth may persist essentially indefinitely since each deposited atom may always 
be able to seek out a desired epitaxial site for incorporation. In some sense the multiple hit 
noise reduction parameter m is equivalent to the dimensionless diffusion length parameter 
l/L because large values of both tend to enhance the layer by layer growth regime. In 
section III, where we present our simulation results, we will see the precise nature of this 
correspondence between m and l/L in our two methods of obtaining layer by layer epitaxy 
in the DT and the WV model. 

The central quantity of interest in layer by layer epitaxial growth is the characteristic 
time t c at which layer by layer growth dies out, i.e. for the deposited average film thickness 
larger than (t c ), measured in lattice units or in monolayers, there are no discernible layer by 
layer growth oscillations. It has been found in earlier numerical simulations of layer by layer 
growth in a variety of contexts that t c obeys an approximate scaling relation with the coarse 
graining parameter m (or l/L as the case may be), and we are interested in investigating 
whether the following scaling relations hold in the limited mobility growth models 



If such scaling relations do hold in our simulations we are interested in obtaining the rela- 
tionship, if any, between the exponents 5 and \i. 

It is in fact fairly straightforward to obtain a relationship between the noise reduction 
parameter m and the surface diffusion length parameter l/L using only dimensional argu- 
ments. In particular, we note that for a <i'-dimensional substrate (d! = 2 for real surfaces 
and d! — 1 for the 1+1 -dimensional growth) one obtains a simple relationship between \i 
and 5 by noting that a surface diffusion length / corresponds tom~ l d ' since there are l d ' 
available sites for a particular deposited particle to incorporate at. This immediately leads 




5 for the long diffusion length method. 



for the noise reduction method 



(1) 



to 
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5 = fid'. (2) 

Within the limited accuracy of our growth simulations (see section III for the results), we 
find the scaling relation defined by Eq. (0) to be valid. 

One secondary theoretical goal of our work is to investigate the extent to which layer 
by layer growth oscillations obey scaling with respect to growth time (or equivalently, the 
average film thickness). Using Eq. (|]) as the theoretical ansatz one can ask whether the 
surface roughness W, defined as the ensemble averaged (over many growth simulations) root 
mean square fluctuation in the interface height, which is oscillatory (and W < 1 implying 
little roughness) in the layer by layer growth regime is a general scaling function of the layer 
by layer growth parameter mor I/I through a dependence of the form : 

W(t) ~ f m (t/m») or h{t/{l/L) 5 ). (3) 

If a scaling form such as Eq. (|3|) holds in the layer by layer growth regime, then Eq. ([I]) for 
t c trivially follows from it — t c being the value of time where layer by layer oscillations cease 
to exist. We could go further in our scaling analyses and ask whether the scaling defined 
by Eq. ([3]) continues to hold (perhaps approximately) well beyond (t > t c ) the layer by 
layer growth regime establishing an approximate scaling relationship between the layer by 
layer growth regime (t < t c ) and the kinetically rough growth regime (t > t c ). Our results 
presented in section III indicate that such an approximate scaling relation does indeed exist 
between the layer by layer growth regime and the kinetically rough growth regime. 

Finally, we note that there have been recent attempts fl32|j33| at developing a theory for 
layer by layer growth oscillations starting from continuum growth equations underlying the 
coarse grained long wavelength behavior of MBE growth. A simple dimensional argument 
32]] , later followed up |33j by a renormalization group approach, leads to the following results 



for MBE growth 



4/3 for d' = 1 

A«H (4) 
2 for d' = 2. 
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The exponents denned by Eq. (|J) correspond to the so-called JT|-§[ conserved fourth order 
nonlinear growth equation or the nonlinear MBE growth equation (which is sometimes also 
referred to as the conserved KPZ equation with nonconserved noise or the Lai-Das Sarma- 
Villain |7|||, LDV, equation) which is given by 

|_ = _ z/4 V 4 / i + A 22 V 2 (V/ i ) 2 + r ? , (5) 

where /i(x, t) is the dynamical height fluctuation variable relative to the average interface (h) 
at the substrate site x (with x is the lateral substrate coordinate), V = is the gradient 
operator along the surface, 7] is the deposition shot noise (which causes the kinetic surface 
roughening), and ^4, A22 are coefficients which in general depend on surface diffusion rate, 
deposition rate, etc. Since the continuum description of the DT and the WV model are 



actually quite complex p3| , |31| , |3^ -p^] , and are likely to be different in different dimensions 
||37|| , it is by no means clear that the exponents (given in Eq. @) corresponding to Eq. 
(i) a Pply without any qualifications to the DT and the WV model (as will be discussed in 
section III where we present our numerical simulations). We therefore also provide below the 
exponent fi for the linear second order Edwards- Wilkinson (EW) growth equation [RSI] which 



applies to the limited mobility Family (F) growth model and may also have significant 



relevance to the DT and the WV model [plTj37. 40 : 



2 for d' = 1 

» = { (6) 
00 for d' = 2. 

The EW equation, whose layer by layer growth exponent in 1+1 and 2+1 dimensions is 
given in Eq. @, is the following : 

— = V2 v 2 h + n. (7) 

We note here that both sets of exponent values given by Eq. (|J) and (|6|) will be relevant in 
our discussion of our simulation results to be presented in the next section. 
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III. RESULTS AND DISCUSSIONS 



We now present our 1+1 (i.e. d! = 1) and 2+1 (i.e. d' = 2) dimensional layer by layer 
growth simulation results for the discrete limited mobility DT and WV models in Figs. 0-|6|. 
In Figs. we present 1+1 dimensional simulation results whereas Figs. |5| and [] give 2+1 
dimensional simulational results for the two growth models. In each figure (to be described 
below) the panel (a) gives the simulated W(t) as a function of growth time t for various 
values of the noise reduction parameter m or the surface diffusion length I with the layer 
by layer oscillations manifestly obvious for larger values of m and (The original DT and 
WV models correspond to the simulation with m = 1 and 1=1, which has no layer by layer 
oscillations by construction.) In panel (b) of each figure we demonstrate our best computed 
scaling collapse of the W — t plots (shown in panel (a)) for various values of m or I with a 
suitable scaling of time t to t/m M or t/(l/L) s as the case may be. For each scaling collapse 
in (b) we try various different values of the exponent fi or 5 to obtain the best statistical 
scaling in the simulated data. The system size used in each simulation is indicated in the 
corresponding figures and the captions. Here we should emphasize that the results shown 
in this paper represent only a typical fraction of our extensive DT and WV layer by layer 
growth simulations. The representative results presented here are of course in complete 
agreement with the full set of our simulation data, and our conclusion is based on a very 
large set of simulation results and not just on the results presented in this paper. 

In Fig. 0we show our 1+1 dimensional DT layer by layer growth simulation results for 
finite surface diffusion length I {m = 1, I > 1) for I = 10, 20, 30, 40, 50. The layer by layer 
oscillations are visually obvious in Fig. 0(a) — the damping time t c increases from t c ~ 10 
for I = 10 to roughly t c ~ 50 for / = 50. In general, the magnitude of the oscillations decays 
exponentially with increasing time, and for t > t c we can only discern the power law increase 
of W ~ where (3 is the growth exponent in the model. In Fig. 0(b) we show our scaling 
collapse of the W{t) data from Fig. 0(a), leading to the exponent value 5 ~ 1.5. Thus, 
t c ~ (//L) 1,5 in d' = 1 DT model. We note that the scaling collapse in the kinetically rough 
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growth regime (i.e. for t > t c ) is not excellent, but it is remarkable that the exponent 5 
which is meaningfully defined only in the layer by layer growth regime continues to provide 
an approximate reasonable description of the kinetically rough growth regime. 

In Fig. ^ we show our 1+1 dimensional DT model layer by layer growth oscillations (Fig. 
0(a)) for the multiple hit noise reduction technique (I = 1, m > 1) for different values of 
the noise reduction factor m. In Fig. |||(b) we show the scaling collapse of the W — t data 
in Fig. [|(a) for various m values. The scaling is excellent with an exponent /i = 1.5. Thus, 
t c ~ m in d! = 1 DT model. 

In Fig. |] we depict our noise reduced layer by layer growth oscillations (Fig. f|(a)) 
and the corresponding scaling collapse of the data for various values of m (with I = 1) in 
the 1+1 dimensional WV model simulations. Again, the scaling exponent /i is found to 
be /x = 1.5 for the best scaling collapse, indicating t c ~ m 1,5 in both DT and WV noise 
reduced models in d! = 1. The finding of the apparent same exponent value /i = 1.5 in 
both DT and WV models in 1+1 dimensional growth is consistent with the fact that the 
effective growth exponent (3 (obtained by plotting In W against In t in the simulated results) 
is almost identical in the two models : From the slope of the log — log plot in Fig. 0(b) we 
obtain (3 ~ 0.338 for the d! = 1 DT model whereas from the slope of the log — log plot in 
Fig. §(b) we obtain (3 ~ 0.339 for the d! = 1 WV model. Thus within the effective time 
and length scales of our simulations the two models (DT and WV) have essentially the same 
effective dynamical universality class, which is consistent with the fact that they have the 
same effective exponent /i = 1.5 in d! = 1. The fact that the asymptotic universality classes 
of the DT and WV models are different even in d! = 1 dimension |H],[37|,f|(J does not seem 
to affect the effective values of \x we obtain in our simulations. 

Before presenting our 2+1 dimensional simulation results in Fig. [5|-[7|we first discuss the 
exponent values 5 and fi, all of which have turned out to be approximately 1.5 in the 1+1 
dimensional DT and WV layer by layer epitaxial simulations. (We do not show here I > 1, 
m = 1 simulation results for the d! = 1 WV model because they are very similar to those 
shown in Figs. 0-^with the same exponent fi ~ 1.5.) First we note that the exponent value 
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is very close to (but somewhat above) the theoretically "expected" exponent \l = 4/3 = 1.33 
predicted in Refs. p2| , |33| assuming that the continuum growth equation for these discrete 



growth models is that given in Eq. (^). We also note that the expected relationship between 
\x and 8, namely 8 = fid' which becomes 8 — fi in d! — 1, is obeyed by our simulation results. 
We also add that the measured exponent fi ~ 1.5 is consistent with other findings in the 
literature [^,^,^|. The cause for our calculated fi (~ 1.5) to be somewhat (by roughly 



10%) higher than the theoretical value (/i = 4/3) is not very clear at this stage. We do not, 
however, believe this discrepancy to be particularly significant because of a number of reasons 
: (1) our scaling collapse are in fact not inconsistent with an exponent of 1.33 although an 
exponent value of 1.5 is definitely a statistically better fit for our scaling collapse; (2) it is 
quite conceivable that there are some systematic finite size and finite time corrections to 
scaling (which are known to be very important in DT and WV models) ; (3) finally, at least 
in the WV model which is definitely known |31|j37|j40[| to asymptotically belong to the EW 



universality class, it is possible that our simulated exponent ft ~ 1.5 is showing some effects 
of the asymptotic universality class since the theoretically expected |3"2"| , |3"B| fi for the linear 
EW equation (our Eq. ([/])) is fi = 2 in d' = 1 dimensions. 

Our d' = 2 dimensional noise reduced (m > 1, / = 1) results for the DT, the WV, 
and the F model are shown in Figs. ||-|7] respectively. These 2+1 dimensional layer by 
layer growth results (as well as the results shown in Fig. 0) in limited mobility models are 
completely new and do not exist anywhere in the literature. We carried out our 2+1 layer by 
layer growth simulations using only the noise reduction technique since our 1+1 dimensional 
results (compare Figs. Q and ^) indicate that the finite diffusion length and the multiple hit 
noise reduction techniques are essentially equivalent and the 2+1 dimensional simulations 
with / > 1 are particularly cumbersome to carry out. Our d' = 2 layer by layer epitaxy 
simulations seem to have produced a few surprising results as discussed below. 

In Fig. ^ we present our d! = 2 layer by layer growth simulations in the DT model 
using the noise reduction technique. The results are depicted in the same manner as in the 
d! = 1 case shown in Fig. § — in particular, Fig. |5|(a) shows the actual layer by layer growth 
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oscillations for various values of m whereas Fig. |5|(b) shows the scaling collapse. It is obvious 
from comparing Fig. and |^ that the layer by layer growth regime is substantially stronger 
in d! — 2 case compared with the d! = 1 case, which is consistent with a much larger value 
of the damping exponent fi ~ 2.5 (compared with 1.5 in d' = 1) in Fig. |5](b) in the d' = 2 
system. Our calculated damping exponent /i ~ 2.5 for the d! = 2 dimensional DT model is 
substabtially (by 25%) higher than the corresponding theoretical prediction |32] , |33|| of /i = 2 
for the LDV equation (Eq.([5|)), which is generally thought to be the continuum description 
for MBE growth. This large discrepancy between our simulated damping exponent ([i ~ 2.5) 
and the theoretical damping exponent (fi = 2) corresponding to Eq. @ in d' — 2 dimensions 
may be a real effect, arising from the recently discovered fact |4]J that the DT model in 



2+1 dimensions has actually a very small (but non-zero) EW V 2 h term in its continuum 
description in contrast to the 1+1 dimensional DT model where the EW V 2 h term strictly 
vanishes (i.e. z/ 2 = in Eq. (|7|)) by virtue of a topological symmetry in the DT model. 
Thus the DT model in 2+1 dimensions may actually asymptotically belong (but with a very 
small value of z/ 2 ) to the EW universality class which according to Eq. ([]) has an infinite 
value of the damping exponent /i. It is therefore possible that the value /i ~ 2.5 in Fig. ^|(b) 
may be indicative of a small correction in the 2+1 dimensional DT model arising from the 
EW term in the continuum equation. More work is needed to conclusively settle this issue. 
This would also explain why the scaling collapse for the damped oscillations in Fig. ||(b), 
particularly for the data in the kinetically rough (t > t c ) growth regime, is not as good as 
the corresponding 1+1 dimensional results shown in Fig. ^| — the asymptotic corrections to 
the LDV equation (Eq. @) arising from a small V 2 /i term (which is present in d' = 2 DT 
simulation results of Fig. [5[ 

The noise reduced d' = 2 WV model simulations, shown in Fig. || are even more 
surprising. The layer by layer growth oscillations in the initial transient time (upto t ~ 10 
or so) are apparent in Fig. ^|(a) for finite values of m although the oscillations are already 
weaker than the corresponding DT results shown in Fig. [5[ This is not expected based on 
the d' = 1 results where Fig. ||| (DT) and Fig. |] (WV) are essentially identical within our 
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simulation sizes and times. Even more surprising is the scaling collapse in the d! = 2 WV 
model shown in Fig. |6|(b) which is clearly not a good scaling behavior — in fact, beyond the 
layer by layer regime (i.e. for t > t c ) there is essentially no scaling behavior in the d! — 2 
WV results. The scaling behavior of the WV model (Fig. 0) in the 2+1 dimensional WV 
model is clearly very different from (and worse than) the corresponding DT results. The 
best damping exponent value for /i we obtain from Fig. |6|(b) is fi ~ 1.5, which is the same as 
the corresponding WV value in dl — 1 as depicted in Fig. £|. We should emphasize that this 
estimate (/i m 1.5) for the damping exponent in the 2+1 dimensional WV model should be 
taken at best as a crude estimate for an effective exponent since there is no scaling behavior 
in the WV data shown in Fig. |6|. Given the very similar d' = 1 behavior in the DT and 
WV models (as shown in Figs. [3] and ^) it is very surprising that the d! = 2 behavior in the 
two models (including the effective damping exponent values fi ~ 2.5 and 1.5 respectively 
for the DT and the WV model) is so completely different. 

What is the explanation for this striking difference in the d! = 2 layer by layer epitaxial 
growth behavior in the DT (Fig. |5|) and the WV (Fig. |B|) model (particularly in view of 
their essentially identical behavior in d' = 1 as seen in Figs. [3] and ^)? The explanation 
actually lies in the recently discovered fact PTlfK] that, while the d! = 1 dimensional WV 



model obeys |[7|,f|(],fyJ] the continuum growth equations given in Eqs. (D and © with 
z/ 4 and A22 7^ and v 2 very small but having a nonvanishing positive value, the d! = 2 
dimensional WV growth model is actually unstable in the sense that the WV morphology 
in the 2+1 dimensional growth forms a regular mounded structure with the mound edges 
having approximately constant slopes. Such an epitaxial mounding instability |42| in the 2+1 



dimensional WV model becomes particularly manifest under the noise reduction technique 
as discussed in details in Ref . jl2] . This unstable mounded morphological growth |52| in the 



noise reduced d! = 2 WV model leads to the peculiar behavior seen for late times in Fig. 
^ where the epitaxial mounding instability prevents the usual layer by layer growth regime 
from behaving in the "usual" manner depicted in Figs. |2|-[5|. 

Finally, in Fig. we show our 2+1 dimensional Family (F) model simulations under noise 
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reduced {m > 1, I = 1) conditions. The F model, which by construction is designed to follow 
exactly the linear second order Edwards- Wilkinson equation (Eq. (|7|)) with z/ 4 = A22 = in 
Eq. ([5]) and 1/2 ^ 0, has very simple growth rules (not shown in Fig. [I]) : Each randomly 
deposited atom seeks to find the site of local height minina as the incorporation site. The F 
model is essentially the discretized version of the EW equation (Eq. (|7|)), and as such has 
the EW dynamical exponents : (3 = and // = 00 in d! — 2. The growth exponent f3 (defines 
the presaturation kinetic roughening of W as W ~ t 13 Jl|-^1) being zero (i.e. W ~ In t) in 
the F model, growth is already very smooth because the roughening is only logarithmic in 
time. In the presence of noise reduction, therefore, the F model (Fig. |7|(a)) shows persistent 
layer by layer growth oscillations in d! = 2 with only logarithmic damping of the oscillation 
induced by kinetic surface roughening. Since the damping exponent /i in the noise reduced 
d! = 2 F model is infinity (Eq. @), we cannot obtain any scaling collapse of the layer by 
layer growth simulation data of Fig. |7|(a) which is obvious in the "scaling plot" shown in 
Fig. 0(b). A very large value of the exponent fi (~ 100, for example) will of course produce 
trivial (and meaningless) data collapse, but we have checked that no finite reasonably small 
(upto ji = 10) value of the damping exponent \i produces scaling in Fg. 0(b). Thus, our 
d' = 2 F results are consistent with the theoretical prediction of fi being infinity in 

the d! = 2 EW equation. We note that we have also carried out d! = 1 dimensional noise 
reduced layer by layer growth simulations (not presented in this paper) in the F model, 
obtaining excellent scaling collapse with the theoretically predicted value of ji = 2. Our 
results for the d' = 1 F model layer by layer growth are consistent with those reported in 



Ref. [43 



IV. CONCLUSION 

We have presented numerical results for extensive computer simulations of various noise 
reduced limited mobility DT 0, WV |TJ, and F |9| growth models in 1+1 and 2+1 
dimensions in order to study the damping of the layer by layer growth epitxy invariably 
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induced by the shot noise inherent in the deposition beam fluctuations. We have used 
both multiple hit noise reduction technique and the long surface diffusion length method 
to obtain the layer by layer growth in our simulations, and have shown that these two 
different techniques for obtaining layer by layer growth are essentially equivalent. Our 
simulation results in general exhibit (with two exceptions noted below) very good scaling 
connecting the layer by layer growth regime with the kinetically rough growth regime. Our 
calculated damping exponents agree well with theoretical predictions where applicable. The 
two exceptions noted above are the 2+1 dimensional WV and the F model where scaling 
fails for different reasons. The 2+1 dimensional noise reduced WV model is known |42[ 



to manifest unstable growth with spectacular epitaxial mounding, which inhibits layer by 
layer growth leading to the failure of scaling collapse. The 2+1 dimensional F model on 
the other hand exhibits very strong and persistent layer by layer growth oscillations (with 
little kinetic roughening) whose damping is expected on theoretical grounds to be extremely 
weak leading to an infinite value of the damping exponent, which is equivalent to saying 
that there is essentially no scaling since in the presence of noise reduction layer by layer 
growth regime lasts forever. 

This work has been supported by the US-ONR. 
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FIGURES 

FIG. 1. Schematic plots showing the diffusion rules for the DT and WV models when the 
diffusion length (1) and noise reduction factor (m) are both one. 

FIG. 2. (a) W-t oscillations for 1+1 DT (L = 1000) with m=l and 1=10, 20, 30, 40, 50 (top 
to bottom), (b) scaling plot of systems in (a) using 5 = 1.5. 

FIG. 3. (a) W-t oscillations for 1+1 DT (L = 1000) with 1=1 and m=l, 5, 8, 10, 15 (top to 
bottom), (b) scaling plot of systems in (a) using \i = 1.5. 

FIG. 4. (a) W-t oscillations for 1+1 WV (L = 1000) with 1=1 and m=l, 5, 8, 10, 15 (top to 
bottom), (b) scaling plot of systems in (a) using fi = 1.5. 

FIG. 5. (a) W-t oscillations for 2+1 DT (L = 10 3 x 10 3 ) with 1=1 and m=l, 5, 8, 10, 15 (top 
to bottom), (b) scaling plot of systems in (a) using jj, = 2.5. 

FIG. 6. (a) W-t oscillations for 2+1 WV (L = 100 x 100) with 1=1 and m=l, 5, 8, 10, 15 (top 
to bottom), (b) scaling plot of systems in (a) using ji = 1.5. 

FIG. 7. (a) W-t oscillations for 2+1 Family (L = 100 x 100) with 1=1 and m=l, 5, 8, 10 (top 
to bottom), (b) scaling plot of systems in (a) using jj, = 0.0. 
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